#!/bin/sh
#
# This is the script for generating files for a specific Dalton test job.
#
# For the .check file ksh or bash is preferred, otherwise use sh
# (and hope it is not the old Bourne shell, which will not work)
#
if [ -x /bin/ksh ]; then
   CHECK_SHELL='#!/bin/ksh'
elif [ -x /bin/bash ]; then
   CHECK_SHELL='#!/bin/bash'
else
   CHECK_SHELL='#!/bin/sh'
fi

#######################################################################
#  TEST DESCRIPTION
#######################################################################
cat > cc2mm_atmvdw.info <<%EOF%
   cc2mm_atmvdw
   ---------------
   Molecule:            Phenol (QM) and 10 water molecules (MM)
   Wave Function:       CC2 / STO-3G
   Molecular Mechanics: SPCE01 model: Point charges and isotropic
                        polarizability introduced in the optimization of 
                        the CC wave function. OLDTG=.TRUE. meaning that
                        point charge relaxed HF orbitals are used.
   Test Purpose:        Run checks energy, dipole moment, quadrupole 
                        moment and ground state excitation energies. Test 
                        the use of atomic van der Wall parameters. 
%EOF%

#######################################################################
#  MOLECULE INPUT
#######################################################################
cat > cc2mm_atmvdw.mol <<%EOF%
ATOMBASIS
QM/MM WATER_1
------------------------
    5    0         1 1.00D-15
        6.0   6    Bas=STO-3G
C          -0.000011        1.204184       -0.241226   0   1
C          -0.000055       -1.179775        1.190184   0   2
C           0.000055       -0.014769       -0.917402   0   3
C          -0.000043        1.225561        1.151838   0   4
C          -0.000082        0.035394        1.873987   0   5
C           0.000048       -1.209960       -0.200235   0   6
        8.0   1    Bas=STO-3G
O           0.000001       -0.102091       -2.281945   0   7
        1.0   6    Bas=STO-3G
H          -0.000071        2.175812        1.667937   0   8
H          -0.000133        0.053625        2.954254   0   9
H           0.000124       -2.144132       -0.743618   0  10
H          -0.000068       -2.110990        1.739952   0  11
H           0.001282        0.792008       -2.639041   0  12
H          -0.000072        2.131743       -0.801250   0  13
   -0.669    10     Bas=MM
O          -0.418929       -1.660265       -4.781945   1   1
O          -4.002114        0.473829       -4.193504   2   1
O           4.342195       -0.872334       -4.151706   3   1
O           4.814551        2.386398       -3.676682   4   1
O          -0.893202        1.678174       -6.674935   5   1
O           0.326288        3.925564       -5.780771   6   1
O          -2.626564       -5.029064       -4.197251   7   1
O           3.402347        1.075004       -6.083000   8   1
O          -2.901438       -0.658868       -6.792613   9   1
O           0.000011       -5.991097       -4.781945  10   1
    0.3345   20     Bas=MM
H          -0.179685       -2.505377       -5.162322   1   2
H          -0.336975       -1.039749       -5.506087   1   3
H          -3.838459        1.406843       -4.056214   2   2
H          -4.906396        0.432394       -4.504486   2   3
H           4.063907       -0.374436       -4.920337   3   2
H           3.934742       -0.420264       -3.412950   3   3
H           5.472248        2.674706       -4.309505   4   2
H           4.254348        1.786929       -4.169589   4   3
H          -0.131046        1.268891       -7.084514   5   2
H          -0.541443        2.450293       -6.231966   5   3
H          -0.072098        3.888013       -4.911281   6   2
H           0.050144        4.770969       -6.134576   6   3
H          -2.610723       -4.593898       -5.049610   7   2
H          -1.858094       -5.599612       -4.205179   7   3
H           2.810562        1.398417       -6.762236   8   2
H           4.194348        0.820678       -6.556489   8   3
H          -2.070827       -0.323058       -7.129474   9   2
H          -3.045368       -0.171928       -5.981257   9   3
H           0.266280       -5.556879       -5.592314  10   2
H           0.658520       -6.672642       -4.647754  10   3
%EOF%
#######################################################################
#  QM/MM INTERACTION INPUT
#######################################################################
cat > cc2mm_atmvdw.pot <<%EOF%
**SYSTP
.NUMMTP  # The number of different MM typesr;ie: SPC, SPCEC1, SPCE01 etc
 1
.TYPE
 0       # Allways start with type 0 which is the QM system!!
.MODEL
 SPC
.CHARGS
 13
 -0.2605
 -0.0165
  0.3157
 -0.0366
 -0.1596
 -0.1918
 -0.4399
  0.0792
  0.0902
  0.1052
  0.0715
  0.3036
  0.1397
.SIGEPS
 13
 6.708867  0.0001115   # Atomic parameters in order: sigma_i, epsilon_i
 6.708867  0.0001115   # Interaction is calculated as 
 6.708867  0.0001115   # Evdw_ij = 4*epsilon_ij*[sigma_ij/R^12-sigma_ij/R^6]
 6.708867  0.0001115   # Construction of two-body parameters are given
 6.708867  0.0001115   # as input in DALTON.INP, *QM3 with the keyword
 6.708867  0.0001115   # .ATMVDW. A multiplication and an additive rule
 5.801443  0.0002709   # is implemented for both parameters:
 4.572986  0.0000478   # SIGADD (+), SIGMLT (*), EPSADD (+), EPSMLT  (*)
 4.572986  0.0000478
 4.572986  0.0000478
 4.572986  0.0000478
 0.000000  0.0000000
 4.572986  0.0000478
*******
.TYPE
 1-10
.MODEL
 SPC_E01
.ALPISO
 1
 9.718
.SIGEPS
 3
 5.954092  0.00024224
 0.000000  0.00000000
 0.000000  0.00000000
*******
**END OF
%EOF%
#
#######################################################################
#  DALTON INPUT
#######################################################################
cat > cc2mm_atmvdw.dal <<%EOF%
**DALTON
.RUN WAVEFUNCTION
*QM3
.QM3
.THRDIP
 1.0D-9
.MAXDIP
 180
.ATMVDW  # (Giving the Van der Waals sigma and epsilon parameters atomic specific)
 SIGMLT  # (Sigma_ij= sqrt[sigma_i*sigma_j]: Option 2: .SIGADD s_ij=1/2(s_i+s_j))
 EPSMLT  # (epsilon_ij=sqrt[epsilon_i*epsilon_j]: Option 2: .EPSADD)
**INTEGRALS
.DIPLEN
.NUCPOT
.THETA
.NELFLD
**WAVE FUNCTIONS
.CC
*SCF INP
.THRESH
1.0D-10
*CC INP
.FREEZE
 7 0
.CC2
.THRLEQ
 1.0D-8
.THRENR
 1.0D-8
.MAX IT
 200
.MAXRED
 600
.MXLRV
 200
*CCSLV
.CCMM
.ETOLSL
 1.0D-7
.TTOLSL
 1.0D-7
.LTOLSL
 1.0D-7
.MXSLIT
 200
*CCFOP
.DIPMOM
.QUADRU
.NONREL
*CCEXCI
.NCCEXC
 2
**END OF
%EOF%
#######################################################################
#  CHECK SCRIPT
#######################################################################
echo $CHECK_SHELL >cc2mm_atmvdw.check
cat >>cc2mm_atmvdw.check <<'%EOF%'
log=$1

if [ `uname` = Linux ]; then
   GREP="egrep -a"
else
   GREP="egrep"
fi

if $GREP -q "not implemented for parallel calculations" $log; then
   echo "TEST ENDED AS EXPECTED"
   exit 0
fi

# MM/MM interaction energy compared:
CRIT1=`$GREP "Eelec \= Sum_n,s\[ \(Q_n\*Q_s\)\/\|R_n - R_s\| \]        \| * (\-|\-0)\.02355955" $log | wc -l`
CRIT2=`$GREP "Epol  \= - 1\/2\*Sum_a\[ Pind_a\*E\^site_a \]          \|   * (\-|\-0)\.00570793" $log | wc -l`
CRIT3=`$GREP "Evdw  \= Sum_a\[ A_ma\/\|R_ma\|\^12 - B_ma\/\|R_ma\|\^6 \] \| * ( |0)\.00465966" $log | wc -l`
CRIT4=`$GREP "E\(MM\/MM\) \= Eelec \+ Epol \+ Evdw                  \| * (\-|\-0)\.02460781" $log | wc -l`
CRIT5=`$GREP "Eelec \= Sum_n,s\[ \(Q_n\*Q_s\)\/|R_n - R_s\| \]        \| * ( |0)\.00435354" $log | wc -l`
CRIT6=`$GREP "Epol  \= - 1\/2\*Sum_a\[ Pind_a\*E\^\(QMclassic\)_a \]   \| * (\-|\-0)\.00006190" $log | wc -l`
CRIT7=`$GREP "Evdw  \= Sum_a\[ A_ma\/\|R_ma|\^12 - B_ma\/\|R_ma\|\^6 \] \| * (\-|\-0)\.00092224" $log | wc -l`
CRIT8=`$GREP "E\(\"QM\"\/MM\) \= Eelec \+ Epol \+ Evdw                \| * ( |0)\.00336940" $log | wc -l`
TEST[1]=`expr $CRIT1 \+ $CRIT2 \+ $CRIT3 \+ $CRIT4 \+ $CRIT5 \+ $CRIT6 \+ $CRIT7 \+ $CRIT8`
CTRL[1]=14
ERROR[1]="THE CLASSICAL MM/MM ENERGY NOT CORRECT"

# QM/MM interaction energy compared:
CRIT1=` $GREP "Epol  \= - 1\/2\*Sum_a\[ MYind_a\*E\^site_a \]         \| * (\-|\-0)\.00572354" $log | wc -l`
CRIT2=` $GREP "( |0)\.00364366.. \| * (\-|\-0)\.00005472.. \| * (\-|\-0)\.00092223.. \| * ( |0)\.00266669.." $log | wc -l`
CRIT3=` $GREP "\-302\.11398359.. \| \-302\.11131689.. \| * (\-|\-0)\.00001561.. \|   ( |0)\.00000000.." $log | wc -l`
TEST[2]=`expr $CRIT1 \+ $CRIT2 \+ $CRIT3`
CTRL[2]=3
ERROR[2]="THE QM/MM ENERGY TERMS ARE NOT CORRECT"

# Dipole moment components compared:
CRIT1=` $GREP "x * (\-|\-0)\.00380000 * (\-|\-0)\.00965863" $log | wc -l`
CRIT2=` $GREP "y * ( |0)\.375402(39|40) * ( |0)\.95417766" $log | wc -l`
CRIT3=` $GREP "z * (\-|\-0)\.22168041 * (\-|\-0)\.56345537" $log | wc -l`
TEST[3]=`expr $CRIT1 \+ $CRIT2 \+ $CRIT3`
CTRL[3]=3
ERROR[3]="DIPOLE MOMENT COMPONENTS ARE NOT CORRECT"

# Quadrupole moment components compared:
CRIT1=` $GREP "1 * \-2\.995115.. * ( |0)\.0116340. * ( |0)\.0036497." $log | wc -l`
CRIT2=` $GREP "2 * ( |0)\.0116340. * 1\.800341.. * \-3\.7964021." $log | wc -l`
CRIT3=` $GREP "3 * ( |0)\.0036497. * \-3\.7964021. * 1\.194773.." $log | wc -l`
TEST[4]=`expr $CRIT1 \+ $CRIT2 \+ $CRIT3`
CTRL[4]=4
ERROR[4]="QUADRUPOLE MOMENT COMPONENTS ARE NOT CORRECT"

# Lowest excitation energies compared:
CRIT1=` $GREP "\^1A   \|    1   \|     ( |0)\.21734..  \|       5\.9142.  \|     47701\.95." $log | wc -l`
CRIT2=` $GREP "\^1A   \|    2   \|     ( |0)\.28834..  \|       7\.8463.  \|     63285\.05." $log | wc -l`
TEST[5]=`expr $CRIT1 \+ $CRIT2`
CTRL[5]=2
ERROR[5]="LOWEST EXCITATION ENERGIES ARE NOT CORRECT"

PASSED=1
for i in 1 2 3 4 5
do 
   if [ ${TEST[i]} -ne ${CTRL[i]} ]; then
     echo "${ERROR[i]} ( test = ${TEST[i]}; control = ${CTRL[i]} ); "
     PASSED=0
   fi
done 

if [ $PASSED -eq 1 ]
then
  echo TEST ENDED PROPERLY
  exit 0
else
  echo THERE IS A PROBLEM 
  exit 1
fi

%EOF%
#######################################################################
